#Removing previous datasets in memory
rm(list = ls())
#Loading the relevant libraries
#install.packages("gghighlight")
library(ggplot2)
library(gridExtra)
library(gghighlight)1 Intro
Open a new Quarto document. Select all text by pressing Ctrl + A (Windows) or Cmd + A (Mac), then press Delete. Next, copy and paste the following preamble.
---
title: "Lab 5"
author: "Your Name"
date: "September 13, 2024"
format:
html:
toc: true
number-sections: true
colorlinks: true
smooth-scroll: true
embed-resources: true
---Render the document and save it under “week5” in your “stats” folder.
2 Background
Before we get into coding, let us go back to a short explanation of how frequncies relate to probability mass functions (pmf) and cummulative distribution functions (cdf).
There are four functions in R for calculating normal distributions: rnorm, pnorm, qnorm, dnorm.
3 Generating Random Values in a Standard Normal Distribution (mean: 0, sd:1)
The rnorm function generates a vector of normally distributed random numbers. We can specify the mean \(\mu\) and the standard deviation \(\sigma\). Remember that that for a standard normal distribution the mean \(\mu = 0\) and the standard deviation \(\sigma = 1\).
Let us see how it works by creating three samples with 10, 100, and 1000 observations.
3.1 Creating a sample of 10 observations
#Set seed allows us to make sure everyone gets the same results
set.seed(123)
#Step1: Creating 10 obs. with mean 0 and sd of 1
generated_data_10obs<-rnorm(10, mean=0, sd=1)
#Step2: Turning the created data into a data frame
generated_data_10obs<-data.frame(generated_data_10obs)
#Step3: Inspecting the first 10 obs.
head(generated_data_10obs, n=10)#Step4: Renaming variable inside the datafarme
names(generated_data_10obs)<-"value"
head(generated_data_10obs, n=10)#Step5: Plotting the data
figure_1<-ggplot(data = generated_data_10obs, aes(x=value))+
geom_histogram()+
theme_bw()+
ggtitle("10 Obs.")
figure_13.2 Creating a sample of 100 observations
set.seed(123)
#Step1: Creating 100 obs. with mean 0 and sd of 1
generated_data_100obs<-rnorm(100, mean=0, sd=1)
#Step2: Turning the created data into a data frame
generated_data_100obs<-data.frame(generated_data_100obs)
#Step3: Inspecting the first 10 obs.
head(generated_data_100obs, n=10)#Step4: Renaming variable inside the datafarme
names(generated_data_100obs)<-"value"
head(generated_data_100obs, n=10)#Step5: Plotting the data
figure_2<-ggplot(data = generated_data_100obs, aes(x=value))+
geom_histogram()+
theme_bw()+
ggtitle("100 Obs.")
figure_23.3 Creating a sample of 1000 observations
set.seed(123)
#Step1: Creating 1000 obs. with mean 0 and sd of 1
generated_data_1000obs<-rnorm(1000, mean=0, sd=1)
#Step2: Turning the created data into a data frame
generated_data_1000obs<-data.frame(generated_data_1000obs)
#Step3: Inspecting the first 10 obs.
head(generated_data_1000obs, n=10)#Step4: Renaming variable inside the datafarme
names(generated_data_1000obs)<-"value"
head(generated_data_1000obs, n=10)#Step5: Plotting the data
figure_3<-ggplot(data = generated_data_1000obs, aes(x=value))+
geom_histogram()+
theme_bw()+
ggtitle("1000 Obs.")
figure_33.4 Comparing the distribution of the three different samples
new_fig<-grid.arrange(figure_1, figure_2, figure_3, ncol=3)set.seed(123)
#Step1: Creating one million obs. with mean 0 and sd of 1
generated_data_milobs<-rnorm(1000000, mean=0, sd=1)
#Step2: Turning the created data into a data frame
generated_data_milobs<-data.frame(generated_data_milobs)
#Step3: Inspecting the first 10 obs.
head(generated_data_milobs, n=10)#Step4: Renaming variable inside the datafarme
names(generated_data_milobs)<-"value"
head(generated_data_milobs, n=10)#Step5: Plotting the data
figure_4<-ggplot(data = generated_data_milobs, aes(x=value))+
geom_histogram()+
theme_bw()+
ggtitle("One Mil. Obs.")
figure_44 Generating the cumulative distribution function (CDF) for a standard normal distribution
The pnorm function generates a vector for a cumulative distribution function of the normal distribution. We can specify the mean \(\mu\) and the standard deviation \(\sigma\). The pnorm will return the probability that the randomly generated variable X takes on a value less than or equal to x:
\(F_X(x) = P (X\leq x)\)
4.1 Sequence from -3 to 3
#99% of a normal distribution falls between -3 and 3
list_seq<-seq(-3, 3, by =0.01)
list_seq_df<-data.frame(list_seq)dummy_df<-pnorm(list_seq, mean = 0, sd = 1)
dummy_df<-data.frame(dummy_df)names(dummy_df)<-"y"
dummy_df2<-cbind(dummy_df, list_seq)figure_4<-ggplot()+
geom_step(data=dummy_df2, mapping=aes(x=list_seq, y=y)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence: -3, 3")
figure_44.2 Highlighting probabilities
This will return the probability that the randomly generated variable X takes on a value less than or equal to -1.
fig4b<-ggplot(data=dummy_df2, aes(list_seq, y)) +
geom_area(fill = "green") +
gghighlight(list_seq < -1)+
theme_bw()
fig4bThis will return the probability that the randomly generated variable X takes on a value less than or equal to 2.
fig4c<-ggplot(data=dummy_df2, aes(list_seq, y)) +
geom_area(fill = "green") +
gghighlight(list_seq < 2)+
theme_bw()
fig4c4.3 Sequence from -10 to 10
#99% of a normal distribution falls between -3 and 3
list_seq<-seq(-10, 10, by =0.01)
dummy_df<-pnorm(list_seq, mean = 0, sd = 1)
dummy_df<-data.frame(dummy_df)
names(dummy_df)<-"y"
dummy_df2<-cbind(dummy_df, list_seq)figure_5<-ggplot()+
geom_step(data=dummy_df2, mapping=aes(x=list_seq, y=y)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence: -10, 10")
figure_5Comparing the two CDFs
new_fig<-grid.arrange(figure_4, figure_5, ncol=2)4.4 Changing size of the randomly generated variable X: increments of 1
list_seq<-seq(-3, 3, by =1)
list_seq[1] -3 -2 -1 0 1 2 3
dummy_df<-pnorm(list_seq, mean = 0, sd = 1)
dummy_df<-data.frame(dummy_df)
names(dummy_df)<-"y"
dummy_df2<-cbind(dummy_df, list_seq)figure_6<-ggplot()+
geom_step(data=dummy_df2, mapping=aes(x=list_seq, y=y)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence: -3, 3\n by 1")
figure_64.5 Changing size of the randomly generated variable X: increments of 0.1
list_seq<-seq(-3, 3, by =0.1)
dummy_df<-pnorm(list_seq, mean = 0, sd = 1)
dummy_df<-data.frame(dummy_df)
names(dummy_df)<-"y"
dummy_df2<-cbind(dummy_df, list_seq)figure_7<-ggplot()+
geom_step(data=dummy_df2, mapping=aes(x=list_seq, y=y)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence: -3, 3\n by 0.1")
figure_74.6 Changing size of the randomly generated variable X: increments of 0.01
list_seq<-seq(-3, 3, by =0.01)
dummy_df<-pnorm(list_seq, mean = 0, sd = 1)
dummy_df<-data.frame(dummy_df)
names(dummy_df)<-"y"
dummy_df2<-cbind(dummy_df, list_seq)figure_8<-ggplot()+
geom_step(data=dummy_df2, mapping=aes(x=list_seq, y=y)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence: -3, 3\n by 0.01")
figure_84.7 Visualizing the results
grid.arrange(figure_6, figure_7, figure_8, ncol=3)5 Generating the inverse of cumulative distribution function (CDF) for a standard normal distribution
The qnorm function generates a vector for the inverse cumulative distribution function of the normal distribution with a mean \(\mu\) and a standard deviation \(\sigma\). It returns the value x such that pnorm(x, \(\mu\), \(\sigma\)) = p
\(P(X\leq x) = F_X(x)\)
Let’s use the same data. First let us examine what we created previously.
head(dummy_df2, n = 10)Because qnorm works in the opposite direction, we will use the probabilities (that the randomly generated variable X takes on a value less than or equal to x) to recreate the values in our sequence
list_seq2<-dummy_df2$y
dummy_df_q<-qnorm(list_seq2, mean = 0, sd = 1)
dummy_df_q<-data.frame(dummy_df_q)
names(dummy_df_q)<-"q"
dummy_df_q2<-cbind(dummy_df_q, list_seq2)Let us now inspect the output
head(dummy_df_q2, n=10)This looks very familiar. Let us now compare it to the previous dataset
head(dummy_df2, n=10)head(dummy_df_q2, n=10)6 Generating the probability density function (PDF) for a standard normal distribution
The dnorm function generates a vector for a probability density function of the normal distribution. We can specify the mean \(\mu\) and the standard deviation \(\sigma\).
dnorm(x, \(\mu\), \(\sigma\)) - probability density function of the normal distribution with mean \(\mu\) and standard deviation \(\sigma\)
list_seq<-seq(-3, 3, by =0.01)
dnorm_df<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df<-data.frame(dnorm_df)
names(dnorm_df)<-"value"
dnorm_df2<-cbind(dnorm_df, list_seq)figure_8<-ggplot()+
geom_step(data=dnorm_df, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()
figure_8Let us see how changing some of these parameters may affect the appearance of our bell curve
6.1 Sequence from -3 to 3
list_seq<-seq(-3, 3, by =0.01)
dnorm_df_3<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df_3<-data.frame(dnorm_df_3)
names(dnorm_df_3)<-"value"
dnorm_df_3_b<-cbind(dnorm_df_3, list_seq)figure_8<-ggplot()+
geom_step(data=dnorm_df_3_b, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence from -3 to 3")
figure_86.2 Sequence from -10 to 10
list_seq<-seq(-10, 10, by =0.01)
dnorm_df_10<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df_10<-data.frame(dnorm_df_10)
names(dnorm_df_10)<-"value"
dnorm_df_10_b<-cbind(dnorm_df_10, list_seq)figure_9<-ggplot()+
geom_step(data=dnorm_df_10_b, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Sequence from -10 to 10")
figure_9Comparing them
grid.arrange(figure_8, figure_9, ncol=2)Let us see what happens if we change the increment
list_seq<-seq(-3, 3, by =0.01)
dnorm_df_inc3<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df_inc3<-data.frame(dnorm_df_inc3)
names(dnorm_df_inc3)<-"value"
dnorm_df_inc3_b<-cbind(dnorm_df_inc3, list_seq)
figure_10<-ggplot()+
geom_step(data=dnorm_df_inc3_b, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Increments of 0.01")
figure_10list_seq<-seq(-3, 3, by =0.1)
dnorm_df_inc3<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df_inc3<-data.frame(dnorm_df_inc3)
names(dnorm_df_inc3)<-"value"
dnorm_df_inc3_b<-cbind(dnorm_df_inc3, list_seq)
figure_11<-ggplot()+
geom_step(data=dnorm_df_inc3_b, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Increments of 0.1")
figure_11list_seq<-seq(-3, 3, by =1)
dnorm_df_inc3<-dnorm(list_seq, mean = 0, sd = 1)
dnorm_df_inc3<-data.frame(dnorm_df_inc3)
names(dnorm_df_inc3)<-"value"
dnorm_df_inc3_b<-cbind(dnorm_df_inc3, list_seq)
figure_12<-ggplot()+
geom_step(data=dnorm_df_inc3_b, mapping=aes(y=value, x=list_seq)) +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Increments of 1")
figure_12Comparing them
grid.arrange(figure_10, figure_11, figure_12, ncol=3)Let’s do some problems now.
7 Normal vs. Non-Normal Distributions
As you might suspect, life expectancy in our dataset is not normally distributed.
#Setting path
setwd("/Users/bgpopescu/Dropbox/john_cabot/teaching/stats/week5/lab/")
#Step1: Loading the data
life_expectancy_df <- read.csv(file = './data/life-expectancy.csv')
library(dplyr)
#Step2: Calculating the mean
life_expectancy_df2<-life_expectancy_df%>%
dplyr::group_by(Code)%>%
dplyr::summarize(life_exp_mean=mean(Life.expectancy.at.birth..historical.))
#Step3: Getting rid of countries
weird_labels <- c("OWID_KOS", "OWID_WRL", "")
clean_life_expectancy_df<-subset(life_expectancy_df2, !(Code %in% weird_labels))
#Step4: Calculating the mean and SDin the dataset
mean_empirical<-mean(clean_life_expectancy_df$life_exp_mean, na.rm = T)
mean_empirical_val<-as.character(round(mean_empirical, 0))
sd_empirical<-sd(clean_life_expectancy_df$life_exp_mean, na.rm = T)
sd_empirical_val<-as.character(round(sd(clean_life_expectancy_df$life_exp_mean, na.rm = T), 2))
#Step5: Creating annotation text
lines<- paste("Mean (\u03BC)", ": ", mean_empirical_val, "\n", "SD(\u03C3)", ": ", sd_empirical_val)
#Step6: Making the histogram for the dataset
y_coord<-17
fig_hist1<-ggplot(data = clean_life_expectancy_df, aes(x=life_exp_mean))+
geom_histogram()+
#Multiplying by 1.5. Otherwise, density will not be visible
geom_density(aes(y=..count..*1.5))+
theme_bw()+
geom_vline(xintercept=mean_empirical, linetype='dashed', col = 'red')+
ggtitle("Life Expectacy as is\n in the Dataset")+
annotate(geom="text",
x=mean_empirical-20,
y=y_coord,
label=lines,
color="red")
#Step7: Generating 235 observations following a normal distribution with same mean and standard deviation
generated_data<-rnorm(235, mean=mean_empirical, sd=sd_empirical)
generated_data<-data.frame(generated_data)
names(generated_data)<-"life_exp_mean"
#Step8: Making the histogram for the generated dataset
fig_hist2<-ggplot(data = generated_data, aes(x=life_exp_mean))+
geom_histogram()+
#Multiplying by 1.5. Otherwise, density will not be visible
geom_density(aes(y=..count..*1.5))+
theme_bw()+
geom_vline(xintercept=mean_empirical, linetype='dashed', col = 'red')+
ggtitle("Life Expectacy generated\n on a Normal Distribution")+
annotate(geom="text",
x=mean_empirical-20,
y=y_coord,
label=lines,
color="red")
#Step9: Comparing the two histograms
grid.arrange(fig_hist1, fig_hist2, ncol=2)How do we do we obtain these graphs? Let us work on the first graph.
7.1 Working with the empirical (actual) data
Let us go back to the life expectancy dataset
Copy and paste the life expectancy dataset from Lab 4 or download it again. You can download the life expectancy dataset below:
7.2 Load the dataset
#Setting path
setwd("/Users/bgpopescu/Dropbox/john_cabot/teaching/stats/week5/lab/")
life_exp_mean <- read.csv(file = './data/life-expectancy.csv')7.3 Calculating the mean
life_exp_mean2<-life_exp_mean%>%
dplyr::group_by(Code, Entity)%>%
dplyr::summarize(life_exp_mean=mean(Life.expectancy.at.birth..historical.))7.4 Clearning the Data
weird_labels <- c("OWID_KOS", "OWID_WRL", "")
clean_life_exp_mean_df<-subset(life_exp_mean2, !(Code %in% weird_labels))7.5 Calculating the mean and SD in the dataset
mean_empirical<-mean(clean_life_exp_mean_df$life_exp_mean, na.rm = T)
mean_empirical_val<-as.character(round(mean_empirical, 0))
sd_empirical<-sd(clean_life_exp_mean_df$life_exp_mean, na.rm = T)
sd_empirical_val<-as.character(round(sd(clean_life_exp_mean_df$life_exp_mean, na.rm = T), 2))7.6 Creating annotation text
lines<- paste("Mean (\u03BC)", ": ", mean_empirical_val, "\n", "SD(\u03C3)", ": ", sd_empirical_val)7.7 Making the histogram for the dataset
y_coord<-17
fig_hist1<-ggplot(data = clean_life_exp_mean_df, aes(x=life_exp_mean))+
geom_histogram()+
#Multiplying by 1.5. Otherwise, density will not be visible
geom_density(aes(y=..count..*1.5))+
theme_bw()+
geom_vline(xintercept=mean_empirical, linetype='dashed', col = 'red')+
ggtitle("Life Expectacy as is\n in the Dataset")+
annotate(geom="text",
x=mean_empirical-10,
y=y_coord,
label=lines,
color="red")
fig_hist18 Working with the made-up (generated) data
8.1 Generating 235 observations following a normal distribution with same mean and standard deviation[235]
generated_data<-rnorm(235, mean=mean_empirical, sd=sd_empirical)
#We can also do:
generated_data<-rnorm(nrow(clean_life_exp_mean_df), mean=mean_empirical, sd=sd_empirical)
generated_data<-data.frame(generated_data)
names(generated_data)<-"life_exp_mean"8.2 Making the histogram for the generated dataset
fig_hist2<-ggplot(data = generated_data, aes(x=life_exp_mean))+
geom_histogram()+
#Multiplying by 1.5. Otherwise, density will not be visible
geom_density(aes(y=..count..*1.5))+
theme_bw()+
geom_vline(xintercept=mean_empirical, linetype='dashed', col = 'red')+
ggtitle("Life Expectacy generated\n on a Normal Distribution")+
annotate(geom="text",
x=mean_empirical-10,
y=y_coord,
label=lines,
color="red")9 Comparing the two histograms
grid.arrange(fig_hist1, fig_hist2, ncol=2)9.1 Generating more observations following a normal distribution with same mean and standard deviation
generated_data_1000obs<-rnorm(1000, mean=mean_empirical, sd=sd_empirical)
generated_data_1000obs<-data.frame(generated_data_1000obs)
names(generated_data_1000obs)<-"life_exp_mean"9.2 Making the histogram for the generated dataset (1000 obs.)
fig_hist3<-ggplot(data = generated_data_1000obs, aes(x=life_exp_mean))+
geom_histogram()+
#Multiplying by 1.5. Otherwise, density will not be visible
geom_density(aes(y=..count..*1.5))+
theme_bw()+
geom_vline(xintercept=mean_empirical, linetype='dashed', col = 'red')+
ggtitle("Life Expectacy generated\n on a Normal Distribution\n 1000 Obs.")+
annotate(geom="text",
x=mean_empirical-10,
y=y_coord,
label=lines,
color="red")
fig_hist3So, we can see that by having more observations, our dataset approaches something that looks like a normal distribution. Note that our mean \(\mu\) and standard deviation \(\sigma\) are the same as in the original data.
10 Problem
Let’s assume that the data in our original empirical dataset is normally distributed with the mean of 61.93 years and a standard deviation of 8.69 years.
10.1 Life expectancy of less than 50 years
What is the probability that a randomly selected country will have a life expectancy of less than 50 years?
ctry_50<-subset(life_exp_mean2, life_exp_mean<50)
ctry_50nrow(ctry_50)/nrow(life_exp_mean2)[1] 0.1171875
The probability that a randomly selected country will have a life expectancy of less than 50 years is 0.12.
Let us also visualize this probability. We will proceed by calculating the frequency table.
#Step1: Rounding the values
new_data2<-life_exp_mean2
new_data2$life_exp_mean_rounded<-round(new_data2$life_exp_mean, 0)
head(new_data2, n=5)#Step2: Creating a frequency table
freq_table<-table(new_data2$life_exp_mean_rounded)
freq_table
38 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
1 2 3 4 4 5 5 3 4 5 12 6 5 6 5 6 7 8 11 10 7 8 14 13 5 15
68 69 70 71 72 73 74 75 77
10 15 14 14 8 6 3 8 4
#Step3: Turning the table into a dataframe
freq_table<-data.frame(freq_table)#Step4: Providing more intuitive names
names(freq_table)[1]<-"life_exp_mean_rounded"
names(freq_table)[2]<-"frequency"
names(freq_table)[1] "life_exp_mean_rounded" "frequency"
#Step5: Turning factor variables into numeric
freq_table$life_exp_mean_rounded<-as.numeric(as.character(freq_table$life_exp_mean_rounded))
str(freq_table)'data.frame': 35 obs. of 2 variables:
$ life_exp_mean_rounded: num 38 43 44 45 46 47 48 49 50 51 ...
$ frequency : int 1 2 3 4 4 5 5 3 4 5 ...
#Step6: Calculaing the PMFs
freq_table$sum_frequency<-sum(freq_table$frequency)
freq_table$pmf<-freq_table$frequency/freq_table$sum_frequency
head(freq_table, n=5)The next step is to create the CDF. A CDF tells us probability that the randomly generated variable X takes on a value less than or equal to x. We will be using the function that we created previously.
#Step7: Defining the function
EvalCdf <- function(t, x) {
count = 0.0
for (value in t){
if (value <= x){
count <- count+1
}
}
prob = count / length(t)
return(prob)
}#Step8: Making list of all rounded ages
list_all_ages_rounded<-new_data2$life_exp_mean_roundedEvaluating the CDF and saving it to a list
#Step9: Evaluating the CDF
xsx<-cbind(sapply(new_data2$life_exp_mean_rounded, function(x){
EvalCdf(list_all_ages_rounded, x)
}))Adding the list to the dataframe:
#Step10: Adding the list to the dataframe
new_data2$cdf <- xsx[,1]Plotting the CDS
#Step11: Plotting the CDS
figure_3<-ggplot()+
geom_step(data=new_data2, mapping=aes(x=life_exp_mean_rounded, y=cdf)) +
geom_bar(stat="identity")+
theme_bw()
figure_3Merging the PMFs to the dataset
#Step12: Selecting only the relevant variables
freq_table_x<-subset(freq_table, select = c(life_exp_mean_rounded, pmf))
head(freq_table_x, n=10)#Step13: Performing the left join
new_data3<-left_join(new_data2, freq_table_x, by=c("life_exp_mean_rounded" = "life_exp_mean_rounded"))#Probability Density
fig4a<-ggplot(data=new_data3, mapping=aes(x=life_exp_mean_rounded, y=pmf)) +
geom_area(fill = "green") +
gghighlight(new_data2$life_exp_mean_rounded < 50)+
theme_bw()+
ggtitle("Probability Density")
fig4a#Cumulative Probability
fig4b<-ggplot()+
geom_step(data=new_data2, mapping=aes(x=life_exp_mean_rounded, y=cdf)) +
geom_area(fill = "green") +
geom_bar(stat="identity")+
theme_bw()+
ggtitle("Cumulative Probability")
#Cumulative Probability < 50
fig4c<-ggplot(data=new_data2, aes(x=life_exp_mean_rounded, y=cdf)) +
geom_area(fill = "green") +
gghighlight(new_data2$life_exp_mean_rounded < 50)+
theme_bw()+
ggtitle("Cumulative Probability")
grid.arrange(fig4a, fig4b, fig4c, ncol=3)If our distribution was normal, it would look different:
It would look something like the following:
If we assume a normal distribution, we can simply get this answer this way:
res<-pnorm(50, mean = mean_empirical, sd = sd_empirical)
res[1] 0.08475828
The probability of getting a country with a life expectancy of less than 50 is 0.12, if we assume a normal distribution.
10.2 Life expectancy of more than 70 years
What is the probability that a randomly selected country will have a life expectancy of more than 70 years?
So remember that the CDF is the probability that the randomly generated variable X takes on a value LESS than or equal to x. We will first subset the data to less than 70 and take the complement.
prob_more_70<-subset(new_data2, life_exp_mean_rounded<70)
1-max(prob_more_70$cdf)[1] 0.2226562
Therefore, probability of getting a country with a life expectancy of greater than 70 is 0.22.
If we assume a normal distribution, we can simply get a similar answer this way:
res<-pnorm(70, mean = mean_empirical, sd = sd_empirical)
1-res[1] 0.1765817
Here is what it looks like
fig4b<-ggplot()+
geom_step(data=new_data2, mapping=aes(x=life_exp_mean_rounded, y=cdf)) +
geom_area(fill = "green") +
geom_bar(stat="identity")+
theme_bw()
fig4c<-ggplot(data=new_data2, aes(x=life_exp_mean_rounded, y=cdf)) +
geom_area(fill = "green") +
gghighlight(new_data2$life_exp_mean_rounded > 70)+
theme_bw()
grid.arrange(fig4b, fig4c, ncol=2)10.3 Percentile
What is the 90th percentile of life expectancy?
We can obtain this by - ordering the ages and calculate the percentile - subset based on CDF and get the minimum
list_ages<-new_data2$life_exp_mean_rounded
list_ages<-list_ages[order(list_ages)]
quantile(new_data2$life_exp_mean_rounded, c(0.90))90%
72
pct_90<-subset(new_data2, cdf>0.9)
#pct_90
min(pct_90$life_exp_mean_rounded)[1] 72
The 90th percentile of life expectancy is 72.
If we assume a normal distribution, we can also use qnorm. This generates a vector for the inverse cumulative distribution function of the normal distribution with a mean \(\mu\) and a standard deviation \(\sigma\). It is not exactly the same, but it is close.
res<-qnorm(0.9, mean = mean_empirical, sd = sd_empirical)
res[1] 73.06729